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An adjoint based technique is applied to a Shallow Water Model in order to estimate 
influence of the model's parameters on the solution. Among parameters the bottom topog- 
raphy, initial conditions, boundary conditions on rigid boundaries, viscosity coefficients and 
the amplitude of the wind stress tension are considered. Their influence is analyzed from 
different points of view. 

Two configurations have been analyzed: an academic case of the model in a square box 
and a more realistic case simulating Black Sea currents. It is shown in both experiments 
that the boundary conditions near a rigid boundary influence the most the solution. This 
fact points out the necessity to identify optimal boundary approximation during a model 
development. 
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ary conditions; Shallow water model. 

1 Introduction 

Thirty years ago model and data were considered as independent one on another. 
Observational data were interpolated on the model grid in order to provide the model 
with the initial conditions, forcings and all the other necessary parameters. However, 
since the pioneering work [1] of Edward Lorenz, we know that a geophysical fluid is 
extremely sensitive to initial conditions. A perturbation of initial state may grow 
exponentially in time limiting the validity of the forecast. This discovery leads to 
understanding that observational data can not be considered as independent of the 
model. We must perform a joint analysis of the model and data in order to choose 
the optimal initial point for the model. 

This become possible by using variational data assimilation technique, first pro- 
posed in [2], [3], which is based on the optimal control methods [4J and perturbations 
theory [5]. This technique allows us to retrieve an optimal data for a given model 
from heterogeneous observational fields ensuring a better forecast. 

However, even now, all other forcings and parameters of the model are obtained 
from data by more or less sophisticated interpolation and they can not be considered 
as optimal for a given model. In the same time, we may suppose, that their influence 
on the models solution is as strong as the influence of initial state. In this case, 
we should also analyze the possibility and utility to apply the data assimilation 
techniques to identify optimal values for all these parameters in order to improve 
the forecast. 

The purpose of this paper is to analyze the sensitivity of a Shallow- Water model 
and, in particular, compare the influence of initial conditions with the influence of 
other parameters. Among these parameters, we consider the boundary conditions on 
the rigid boundaries, bottom topography, empirical coefficients like reduced gravity, 
forcing amplitude and dissipation. 
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2 Sensitivity and control of the boundary conditions 



We shall focus our attention on the boundary conditions because (as we shall see 
later) they represent the most unusual control variable. 

However, as it has been noted in [6j, particular attention must be paid to the 
discretization process which must respect several rules because it is the discretiza- 
tion of the model's operators that takes into account the set of boundary conditions 
and introduces them into the model. Consequently, instead of controlling bound- 
ary conditions them-self, it may be more useful to identify optimal discretization of 
differential operators in points adjacent to boundaries because this is more general 
case. Indeed, boundary conditions participate in discretized operators, but consid- 
ering the discretization itself, we take into account additional parameters like the 
position of the boundary, lack of resolution of the grid, etc. 

Boundary conditions are usually introduced into the model by a particular 

discretization of operators near the boundary. For example, taking into account 

the condition uq = we can calculate the derivative at the point x = h/2 as 
du _ ui 
^ 1/2 "TT' 

In this paper, we shall write the approximation of the derivative in a general 
form 



du 
dx 



1/2 ^ 



Coefficients ao and ai will be used as controls. That means we shall let them vary 
in the data assimilation procedure in order to find an optimal pair that realizes the 
minimum of the cost function. 

2.1 Example: one-dimensional wave equation 

In order to understand what happens when the data are assimilated to control 
the boundary conditions, we propose to take a look on a scholar example: one- 
dimensional wave equation written for u = u{x,t) and p = p{x,t) in the following 
way: 

du dp ^ 

dt dx 

I - « 

This equation is defined on the interval < x < 1 with boundary conditions pre- 
scribed for u only: 

u{0,t) =u{l,t) =0 (2) 
Initial conditions are prescribed for both u and p 

u{x, 0) = u, p{x, 0) = p (3) 

The equation is discretized on a regular grid that is somewhat similar to Arakawa 
C grid [7] in two dimensions: 



Pl/2 ^ P3/2 ^ P 5/2 ^ PiV-5/2 PiV-3/2 PiV -1/2 

ui U2 113 "N-3 "N-2 "N-1 
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Discrete derivatives of u and p are defined as follows 

' dp\ _ Pi+l/2 -Pi-l/2 fdu\ _ Uj+i - Uj 

dxJi h ' \dx)i+ii2 h 

in all internal points i.e. 2 < i < iV - 2 for and 1 < i < iV - 2 for 

Near the boundary, at points 1/2, 1, — 1, — 1/2, we write the derivatives in a 
general form, like 

dp\ _Oilpij2 + o?iP-i/2 (du\ _a^ + a>i 



dx J I h ' \dx J h 

considering oq and ai as the control coefficients. Leap-frog scheme was used for 
time stepping. 

We introduce the simplest cost function that represents the distance between 
the model solution and observation at time t: 

T 1 

X{a) = j J u{a,x,t) -u"^%x,t)f + {p{a,x,t) -p"^%x,t)fdxdt (6) 



and we calculate its gradient using the adjoint to the derivative of the solution with 
respect to control coefficients aP,a^: 

J { du{t),p{t) y ( n(a,x,t)-n°''^(x,t) V 

V I da ) [ p(a,x,t)-p°^^(x,t) 
^ ^ 

Once we prescribe the initial conditions for the equation 

u{x,0) = sin^kiTx) p{x,0) = cos{k'iTx), 
we can calculate its exact solution: 

Uexact{x,t) = V2sin{kiTt—TT / 4:)sin{kTTx) , Pexact{x,t) = —V2cos{kTTt—7r/4)cos{kTrx). 

The exact solution is used as artificial observational data in this example. We 
perform the minimization of the cost function The minimization procedure 

used here was developed by Jean Charles Gilbert and Claude Lemarechal, INRIA 
[8]. The procedure uses the limited memory quasi-Newton method. 

The difference between the models solution and the exact one is shown in figUl 

We see that optimal discretization of derivatives near the boundary brings the 
solution much closer to the exact solution, but the set of optimal coefficients a does 
not approximate a derivative: 

du\ _^o^3 ^i-^o (dp\ _ 3.014P3/2 - 2.828pi/2 
dxJi/2 ' /i ' \dx J I h 

Neither expression for nor for ^ has any reasonable order of approxima- 
tion. The first one is of order, the second is of -1 order. Moreover, while we get 
always the same formula for ^ , approximation of the derivative of p varies in dif- 
ferent assimilation experiments. Assimilations performed with different assimilation 
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Figure 1: Difference between the models solution and the exact one: classical BC 
-dashed line, optimal BC - solid line 

windows, for example, result in different coefficients for In fact, any combination 

V 

Oq , a 



in ([5]) may be found as the result of assimilation under condition 

a\ = -1.104ag - 0.107. (9) 



This linear relationship has been obtained experimentally performing assimilations 
with all assimilation windows in range from 600 to 2400 time steps (with the time 
step equal to 1/120 of the time unit). 

To explain this strange result, we analyze the numerical solution of the equation. 
It is well known, the principal numerical error of the scheme is a wrong wave velocity. 

The wave speed, that must be equal to 1, is replaced by /3 = 2Tsin{kh]2) ^'^^^^ 

depends on the time step r and the grid step h. For the given parameters {k = 

3, h = and r = j^U-'' ^^^'^^ ™ wave velocity is equal to 3.09 x 10^'^ 

The data assimilation and control of the boundary derivatives can not modify 
numerical wave velocity. The only way for this control to get a better solution 
consists in modifying the length of the interval. A numerical wave with wrong 
velocity will propagate on the interval with wrong length. But the length of the 
interval is adapted by data assimilation in order to ensure the wave with numerical 
velocity propagates the modified interval in the same time that the exact wave 
propagates the exact interval. So far, the control can not correct the error in the 
wave velocity, it commits another error in length in order to compensate the first 
one as it is illustrated in figi2j 

Pl/2 P3/2 PS/2 
I X 1 X \ X 



"0 



Pl/2 P3/2 PS/2 

1 \ ¥r 



Figure 2: Modification of the intervals length. 

Non uniqueness of optimal g?y and Oq can be explained if we take into account 
that -p has also a form of cosine of Svrx. Hence, at any time = ^(i) cos(37r/i/2) 
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and = ^(t) cos(97r/i/2) with some A depending on time. Their hnear combina- 
tion a'lpi/2 + 00^3/2 can vanish if 

"^""4cos2(A:7rV2)-3- 

Consequently, all couples , Oq belonging to the line that passes by the point 

On = —1.023 ,a? = 1.023 with tangent — -; ^, = —1.108 produce the 

" ^ 4cos^(37r/i/2) - 3 

same derivative. This line coincides withing accuracy of computation with the set ([9]) 
obtained numerically. Any point on this line gives coefficients ap that theoretically 
provide the same value of the derivative and the same value of the cost function. 

Of course, in this simple example we can avoid the ambiguity in the solution: it 
is sufficient to control only one coefficient rather than two. But, in more complex 
problems, it may be difficult to locate and avoid the presence of kernels. 

Consequently we can say that the data assimilation allows to place the boundary 
in the optimal position resulting in a solution closer to the exact one. Boundary con- 
trol allows to compensate numerical errors committed in the interior of the domain, 
but it may be difficult to understand the physical meaning of optimal coefficients a 
and non-null kernels may exist leading to non unique result. 

More details of this study can be found in [9] 

2.2 Shallow Water Model 

In this paper we consider a shallow-water model written in a conservative form: 
dhu ^, d f , 2 ■ T /, , du\ 



fhv-—[ hu^ + gh{h - H) - fih- 



dt dx\ dx J 



huv - nh^^ - ahu + TQT^, 
oy\ oyJ 

' t ■ 111 w r \ I 1 



^ /iv + gh{h - H) - f^^'Q:j^j ~ '^^^ + '''oTy, (11) 



dh dhu dhv 



dt dx dy 

where hu{x, y, t) and hv{x, y, t) are two flux components that represent the product 
of the velocity by the ocean depth, h{x,y,t), that corresponds to the distance from 
the sea surface to the bottom of the ocean. The sea surface elevation is represented 
by the difference h{x,y,t) — H{x,y), where H{x,y) is the bottom topography. The 
model is driven by the surface wind stress with components Tx{x, y, t) and Ty(x, y, t) 
normalized by tq and subjected to the bottom drag that is parameterized by linear 
terms ahu and ahv. Horizontal eddy diffusion is represented by harmonic operators 
div{fj,h\/u) and div{fih'Vv). Coriolis parameter is represented by the variable f{y) 
that is equal to fo+Py assuming /3-plane approximation. Parameter g is the reduced 
gravity. The system is defined in some domain with characteristic size L requiring 
that both hu and hv vanish on the whole boundary of 0,. No boundary conditions 
is prescribed for h. Initial conditions are defined for all variables: hu, hv and h. 

As usual, initial conditions are considered as the control parameter of the model 
in this paper. We study the sensitivity of the model to its initial point and as- 
similate data to find its optimal value. However, in addition to initial conditions, 
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all other parameters of the model, and namely the discretization of operators near 
the boundary, its bottom topography H{x,y), scalar coefficients fJ,,cr,g and tq, are 
also considered as control variables. All of them are allowed to vary in the data 
assimilation procedure in order to bring them to their optimal values. 

We discretize all variables of this equation on the regular Arakawa's C-grid |10] 
with constant grid step 6x = ^ in both x and y directions. Discretizing the system 
(llip . we replace the derivatives by their finite difference representations and Dy 
and introduce two interpolations in x and y coordinates Sx and Sy. Interpolations 
are necessary on the staggered grid to calculate the variable's values in nodes where 
other variables are defined. The discretized system (fTTj) writes 



^^"^ fSxSyhv + Dx ( ^-^4^ + ghih - Hix, y)) - f^hD,^ ) + 



({SyhuSxhv) (Q Q u\r> h 1 

^ S S h f^i^^^y'^J^y^ J = "^^^ + '^oTx, 

— + fSySxhu + Dx^ ^^^^^ - KSxSyh)Dx—) + 

+ Dy (^i^y!^ + gh{h - H{x, y)) - l^hDypj-^ = -ahv + ToTy, (12) 



Drhu — D.,hv. 



dh 
'dt 

Discretized operators Dx, Dy and S^, Sy are defined in a classical way at all internal 
points of the domain. For example, the second order derivative and the interpolation 
operator of the variable hu defined at corresponding points write 

{Dxhu)i_y2,j-i/2 = for z = 2, . . . , - 1, 

h \ huij_i/2 + hUi_ij_i/2 p . „ AT ^ 
{Sxhu)i_i/2j-i/2 = ^ for ? = 2, ...,A^-1. (13) 

Discretization of operators in the directly adjacent to the boundary nodes are 
different from (jl3p and represent the control variables in this study. In order to 
obtain their optimal values assimilating external data, we suppose nothing about 
derivatives and interpolations near the boundary and write them in a general form 

{Dxhu)i/2j-l/2 = Oo + 

, . gfe^ , ai" ^no,j-i/2 + «2" huij_i/2 
[Sxhu)i/2j-i/2 = Oq + 2 ^ ^ 

This formula represents a linear combination of values of hu at two points ad- 
jacent to the boundary with coefficients a. The constant ao may be added in some 
cases to simulate non-uniform boundary conditions like hu{0, y) = ^ 0. 

We distinguish a for different variables and different operators allowing different 
controls of derivatives because of the different nature of these variables and differ- 
ent boundary conditions prescribed for them. It is obvious, for example, that the 
approximation of the derivative Dx in the first equation may differ from the approx- 
imation of Dx in the third one. Although both operators represent a derivative, 
boundary conditions for hu and h are different, these derivatives are defined at dif- 
ferent points, at different distance from the boundary. Consequently, it is reasonable 
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to let them be controlled separately and to assume that their optimal approximation 
may be different with distinct coefficients a^^" and a^^^ . 

Time stepping of this model is performed by the leap-frog scheme. The first 
time step is splitted into two Runge-Kutta stages in order to ensure the second 
order approximation. 

As well as before, the approximation of the derivative introduced by (jl3p and 
(I14p depends on variables a. These variables are added to the set of control variables 
enumerated above. Operators are allowed to change their properties near boundaries 
in order to find the best fit with requirements of the model and data. To assign all 
control variables we shall perform data assimilation procedure and find their optimal 
values. Variational data assimilation is usually performed by minimization of the 
specially introduced cost function. The minimization is achieved using the gradient 
of the cost function that is usually determined by the run of the adjoint to the 
tangent linear model. 

To define the cost function we introduce dimensionless state vector (/> that is 

composed of three variables of the model = {whuhu,Whvhv,Whhy weighted by 

coefficients w. These weights are used to normalize values of the flux components 

1 1 
by Whu = Whv = T7 — / and the Sea surface elevation by = -rj-- The distance 

-nov9-"0 -"0 
between the model solution and observations is defined as the Euclidean norm of 

the difference 

e = e{Hp,t)) = j2('pk-<Pt?= (15) 

k 

i,j i,j i,j 

In this expression, we emphasize implicit dependence of ^ on time and on the set of 
the control parameters p that is composed of 

• the set of initial conditions of the model (po = {hu \t=o, hv \t=o, h \t=o}, 

• the set of the coefficients a that controls the discretizations of operators near 
the boundary, 

• the bottom topography H{x, y) 

• four scalar parameters fj, /U, 5, To . 

Taking into account the results obtained in [Tl] , we define the cost function as 

T 

ap) = J te{<p{p,t))dt (17) 



that gives higher importance to the difference at the end of assimilation interval. 

It should be noted here, that this cost function can only be used in the case of 
assimilation of a perfect artificially generated data. When we assimilate some kind 
of real data that contains errors of measurements and is defined on a different grid, 
we should add some regularization term to the cost function (like the distance from 
the initial guess) and use some more appropriate norm instead of the Euclidean one 
(see, for example for details). 
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The reth component of the gradient of the cost function can be calculated as the 
Gateaux derivative of an implicit function: 







T 

^Jt(j:{^k-€'n^)dt (18) 





because the derivative j^j^ can easily be calculated from p5|) : = 2{(f)k — (t>k) 



The second term in (jl8p , -jv-^ , represents the matrix of the tangent linear model that 
relates the perturbation of the parameter pn and the perturbation of kth component 
of the model state vector cpk- This relationship, of course, is assumed in the linear 
approach, that means it is only valid for infinitesimal perturbations. 

In the classical case, when initial conditions are considered as the only control 

variable, the derivative ^'^^ = classical tangent model that describes 

the temporal evolution of a small error in the initial model state. The matrix is a 
square matrix that is widely studied in numerous sensitivity analyses. Its singular 
values at infinite time limit are related to well known Lyapunov exponents that 
determine the model behavior (chaotic or regular) and the dimension of it's attractor. 

In our case, the matrix ^'^^ is rectangular. It describes the evolution of an 

infinitesimal error in any parameter (including initial state). However, we can study 
it's properties in the similar way as we do with the classical tangent linear model. 
Its structure and composition is described in [11] for the case of using coefficients a 
as control parameters and in [T3] for the case when the bottom topography is used 
to control the model's solution. ^ 

product. To calculate this product directly we would have to evaluate all the ele- 
ments of the matrix. This would require as many tangent model runs as the size 
of the state vector is. So, instead of the tangent model, we shall use the adjoint 
one that allows us to get the result by one run of the model. Backward in time ad- 
joint model integration that starts from [(p — cj)"^^) provides immediately the product 



-^j{(l>- <P"'") which is exactly equal to {(p - (t)"^")^ in ([T8 



Using these notations, we write 



vx=2 / tm}^Y{<p{p,t)-r'mt (19) 



obs I 







\ dp J 



where the expression in the integral is the result of the adjoint model run from t to 
starting from the vector {(t>{p,t) — 4>°''^{t))- 

Tangent and adjoint models have been automatically generated by the Tapenade 
software [14j.[15j developed by the TROPICS team in INRIA. This software analyzes 

the source code of the nonlinear model and produces codes of it's derivative g and 
of the adjoint ( ^ 
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This gradient is used in the minimization procedure that is implemented in order 
to find the minimum of the cost function: 



I{p) = minX(p) (20) 

Coefficients p are considered as coefficients achieving an optimal parameters for the 
model. As it has been already noted, the set of parameters p is composed of the set 
of initial conditions of the model (po, the set of the coefficients a that controls the 
discretization of operators near the boundary, the bottom topography H{x, y) and 
four scalar parameters a,^,g,TQ. We shall minimize the cost function controlling 
either the total set of available parameters p or any possible subset, comparing the 
efficiency of the minimization. 

We use the minimization procedure developed by Jean Charles Gilbert and 
Claude Lemarechal, INRIA [8|. The procedure uses the limited memory quasi- 
Newton method. 

In addition to the data assimilation, we perform also the sensitivity study of the 
model solution to parameters enumerated above. We are looking for a perturbation 
in the model's parameters 6p that, for a given small norm, maximizes the norm of 
the perturbation of the solution at time t. 



X{t) = max " 7 " (21) 

Sp 



We can note that we already have all the necessary software to estimate \{t). Tan- 
gent linear model { ^'^^ ^ allows us to calculate 5(j){t) = { ^'^^ ^^P- Using the 

scalar product that corresponds to the norm in the definition of the distance ^ (|16|) . 
we can write 

m)\,^ (dm 



X{t) = max — = max ^ — = 

<^op, dp:$> -^op, dp:$> 

^6p, 5p» 

This expression is a well known Rayleigh-Ritz ratio which is equal to the largest 
eigenvalue of the problem 

So far, we need just the maximal eigenvalue and the matrix of the problem is a 
self-adjoint positive definite matrix, we can solve the problem (j23p by the power 
method performing successive iterations 

\ op J \ op J 
= ,.,xx * / ' ^0 = random vector 

In the limit, the denominator of the right-hand-side tends to the largest eigenvalue 
and i?„ — to the corresponding eigenvector of the matrix. The principal advantage 
of this method consists in the fact that we do not need to calculate the matrix itself. 
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we just need a matrix-vector product. So far, we have both codes for the tangent 
and adjoint models, we can successively run these models and get the left-hand side 
of (i23]). 

We should note here that when the initial conditions of the model are used as 
the control parameters (i.e. 5p = (5(/)(0)), the sensitivity characteristics A(t) are all 
close to one when t — > 0. It is evident because the perturbation has no time to be 
transformed by the model's dynamics and we get (5</>(t) \t >o= <^0(O) = Sp. 

When any other model parameter is used as the control and the error growing 
time is small, all X{t) are vanishing. This is also clear: the model's dynamics has 
no time to transmit the perturbations from the parameters to the solution. The 
perturbation of the solution remains, consequently, close to zero as well as the value 
of A(t) \t^o= 0. 

In order to make the behavior of the sensitivity characteristics uniform with 
different parameters, we shall use \{t) — 1 every time when the initial model's state 
is considered as the control parameter. 



3 Configurations. 

3.1 Model in a square box. 

We start from the data assimilation in frames of the very well studied "academic" 
configuration. Several experiments have been performed with the model in a square 
box of side length L = 2000 km driven by a steady, zonal wind forcing with a 
classical sinusoidal profile 

2n{y - L/2) 

Tx = To cos 

that leads to the formation of a double gyre circulation [16]. The attractor of the 
model and the bifurcation diagram in a similar configuration has been described 
in [T7]. Following their results, we intentionally chose the model's parameters to 
ensure chaotic behavior. The maximal wind tension on the surface is taken to 

be tq = 0.5 '^^^■f . The coefficient of Eckman dissipation and the lateral friction 

cm 

2 

coefficient are chosen as cr = 5 x lO^^s^^ and ^ = 200^^ respectively. 

As it has been already noted, the Coriolis parameter is a linear function in y 
with /o = 7 X 10~^s~^ and /3 = 2 x 10~^^(ms)~^. The reduced gravity and the 

depth are respectively equal to g = 0.02^, Hq = 1000m. 

s 

The resolution of the model in this section is intentionally chosen to be too coarse 
to resolve the Munk layer [18] that is characterized by the local equilibrium between 
the /3-effect and the lateral dissipation. Its characteristic width is determined by the 

Munk parameter d = 2(^1 which is equal to 42 km in the present case. The 

model's grid is composed of 30 nodes in each direction, that means the grid-step is 
equal to 67 km, that is more than the Munk parameter. Thus, there is only one grid 
node in the layer and the solution exhibits spurious oscillations near the western 
boundary due to unresolved boundary layer. 

Artificial "observational" data are generated by the same model with all the 
same parameters but with 9 times finer resolution (7.6 km grid step). The fine 
resolution model, having 7 nodes in the Munk layer, resolves explicitly the layer and 
must have no spurious oscillations. All nodes of the coarse grid belong to the fine 
grid, consequently, we do not need to interpolate "observational" data to the coarse 
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grid. We just take values in nodes of the high resolution grid that correspond to 
nodes on the coarse grid. 

The model on the fine grid has been spun up from the rest state during 3 years. 
The end of spin up was used as the initial state for the further integration of the 
model. From the result of this integration we have extracted values of all three 
variables at all grid points that belong to the coarse grid (as it has been noted, the 
grids have been chosen so, that all grid points of the coarse grid belong to the fine 
grid). This set is used as artificial observations in the following experiments. 

So far the model is nonlinear with intrinsicly instable solution, there is no hope 
to obtain close solutions in long time model runs because any difference (even in- 
finitesimal) between two models grows exponentially in time. Consequently, we have 
to confine our study to the analysis of a short time evolution of the model's solution 
simulating the forecasting properties of the model. 

As the initial guess for the initial conditions we use the state vector of the 
high resolution model reduced on the coarse grid. This state is also used as the 
initial conditions in all other assimilation experiments with other control parameters. 
Noted above values of the model's parameters (flat bottom topography, linear in y 
Coriolis parameter and scalar parameters (/x, cr, tq, g) are used as the initial guess 
in the experiments that control these parameter, otherwise we simply use these 
parameters in the model. 

3.2 Model of the Black Sea. 

In this section we use the same model, but all the parameters are defined to de- 
scribe the upper layers circulation of the Black sea. Configuration of the model and 
observational data have been kindly provided by Gennady Korotaev from the Ma- 
rine Hydrophysical Institute, National Academy of Sciences of Ukraine, Sevastopol, 
Ukraine. This configuration is described in |19j . 

The model grid counts 141 x 88 nodes that corresponds to the grid box of 
dimension 7860 m and 6950 m in x an y directions respectively. 15 minutes time step 
is used for integration of the model. The Coriolis parameter is equal to /o = 10~^s~^ 
and /3 = 2 X 10~^^ (ms)~^. Horizontal viscosity is taken as /i = 50m^s~^. Using 
a typical density difference between upper and underlying layers of ^.Ikg/m?, and 
unperturbed layer thickness of Hq = 150m, the Rossby radius of deformation is 
estimated at about 22 km and the reduced gravity value g = 0.031m/s^. The grid 
therefore resolves the mesoscale processes reasonably well. 

The model has been forced by the ECMWF wind stress data, available as daily 
averages for the years 1988 through 1999. Dynamical sea level reconstructed in 
[20] was used as observational data in this section. These data have been collected 
in ERS-1 and TOPEX/Poseidon missions and preprocessed by the NASA Ocean 
Altimeter Pathfinder Project, Goddard Space Flight Center. Observational data 
are available from the 1st May 1992 until 1999. These data have been linearly 
interpolated to the model grid. 

So far the sea surface elevation is the only observational variable available in this 
experiment, we put Whu = Whv = in (fTB]) . Consequently, the difi^erence between 
the model's solution and observations is calculated taking into account the variable 
h only. 

As it has been already noted, absence of observational data for the velocity 
fields brings us to modify the cost function. We have to add the background term 
in the cost function in order to require the velocity field to be sufficiently smooth. 
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Otherwise, lack of information about velocity components in observational data 
would result in a spuriously irregular fields obtained in assimilation. To ensure 
necessary regularity of hu and hv we add the distance from the initial guess to 
the cost function ()17p . In order to emphasize the requirement of smoothness, this 
distance is measured as an enstrophy of the difference between the initial guess and 
current state: 



-'Smooth — / , 



dx dy 



(24) 



where hu^ , hv^ denote flux components of the initial guess of the minimization 
procedure. 

Moreover, using real observational data requires to add at least one another 
term to the cost function. One can see in the Figure 2 in [20], spatially averaged sea 
surface elevation of the Black sea exhibits a well distinguished seasonal cycle. That 
means the mass is not constant during a year, it decreases in autumn and increases 
in spring. Consequently, if we assimilate data during a short time (a season or 
less), we assimilate also the information about the mass flux specific for this season. 
This flux can not be corrected later by the model because the discretization of 
operators near the boundary (that controls the mass evolution) is obtained once for 
all seasons. The mass variation of the Black sea reaches 25 centimeters of the sea 
surface elevation. Assimilating data within one season may, consequently, result in 
a persisting increasing or decreasing of the seal level of order of 50 cm per year. To 
avoid this spurious change of the total mass, we must either take the assimilation 
window of at least one year, or prescribe the mass conservation to the model's 
scheme. One year assimilation window is computationally expensive and is not 
justified by the model's physics. On the other hand, prescribed mass conservation 
removes just the sinusoidal seasonal variation, allowing us to keep all other processes 
and to choose any assimilation window we need. 

To correct the mass flux of the model, we add the following term to the cost 
function 

Tmass = J - hji^))) dt (25) 



Similarly to (I24p . this term also ensures the regularity of the solution. It can be 
noted here that other terms may be added to the cost function in order to make a 
numerical scheme energy and/or enstrophy conserving, but we do not use them in 
this paper. 

The total cost function in this section is composed of three parts: (I17p . (|24p and 
(I25|): 

'^total = 2^ + ll^smooth + l2'^mass (26) 

Coefficients 7 are introduced to weight the information that comes from observa- 
tional data (with X) and an a priori knowledge about mass conservation and regu- 
larity of the solution. 

This modification of the cost function results, of course, in additional terms in 
the gradient: 

y^totai = VI+2ji(^D;Dy{hu-huo) + D*Mhv-hvo)^ +272 ^ (^?7ij(t) -ryij(O) 



(27) 
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The model is spun up from the beginning of 1988 to May 1992 using the wind 
tension data on the surface. The state corresponding to the 1st of May 1992 12h 
GMT is used as the initial guess in the data assimilation procedure controlling initial 
conditions of the model. The assimilation controls the initial conditions (po only with 
the assimilation window T = 1 day and the regularization parameter 71 = 0.04. Such 
a short window was chosen in order to get almost instantaneous state of the model 
to be used in further experiment as an initial state. 

The behavior of the model solution is not chaotic in this configuration. Vari- 
ability of the model is generated directly by the variability of the wind stress on 
the surface. Consequently, we can compare particular trajectories of the model on 
any time interval because their evolution is stable without exponential divergence. 
Thus, we can hope that assimilating data in a relatively short window allows us to 
bring the model's solution closer to observation for a long integration period. 

The minimization of the cost function has been accompanied by the mass pre- 
serving correction (j25p with 72 = 0.01. 



4 Sensitivity analysis. 

The flexibility of the model is illustrated in figl3j We perform the data assimila- 
tion experiment in two configurations using parameters described above as initial 
guess. Due to high CPU time of the data assimilation, we limit the number of iter- 
ations of the minimization procedure by 20. Thus, we have similar and reasonable 
computational cost in each experiment. 

In both configurations we examine the evolution of the distance " model-observations' 
^(t) during assimilation and after the end of assimilation. Assimilation window has 
been chosen as 5 days in the square box configuration and T = 30 days for the Black 
sea model. The distance is examined over longer intervals: 20 days in the first case 
and 1 year in the second one. 




Figure 3: Distance between the model solution and observations for the model in 
the square box (left) and the model of the upper layer of Black Sea (right). 

Analyzing the difference between the model solution and observations shown in 
figiSl we see that in the assimilation window the model is almost equally flexible 
with respect to both initial and boundary conditions. Data assimilation allows us 
to reduce the distance between the model solution and observations at the end of 
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the window approximately twice in both configurations. The only difference that 
can be seen in the assimilation window is that non-optimal initial point leads to the 
spurious oscillations of the solution. These oscillations occur in both configurations 
and show us the necessity to identify the optimal initial point. 

However, the influence of parameters is significantly different beyond the win- 
dow. While the solution with optimal initial point tends towards the solution ob- 
tained without any data assimilation, optimal set of boundary conditions ensures 
a new solution that is much closer to observational data. That means the control 
of boundary conditions allows us to improve a long-range forecasting quality of the 
model. 

The third way of the sensitivity analysis consists in solving of the eigenvalue 
problem (I23|) and analyzing A(t) on different scales of error growing time from about 
10 minutes (10"'^ day) to approximately one year. As it has been already noted, 
X{t) — 1 is plotted in the case when initial conditions are considered as the parameter. 
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Figure 4: Finite time Lyapunov exponents on short (above) and long (below) time 
scales: model in the square box (left) and the model of the upper layer of Black Sea 
(right). 

Analyzing the figure figiH we can see that three time scales can be clearly 
distinguished for the sensitivity characteristics of the model in both configurations. 
The first, short time scales, approximately from up to 2-3 hours is characterized by 
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the linear growth of A(t). Indeed, the model behaves as a linear model on this scales, 
the model's solution can be well approximated by just one step of the numerical time 
scheme. 

The second time scale that can be distinguished in the figure figlH corresponds 
to error growing times from 2-3 hours to 10 days. On these time scales we see 
slower growth of the sensitivity characteristics A(t) and, sometimes, no growth at 
all. These time scales are characterized by the modification of the stable-instable 
subspaces of the model. Instable space on short time scale is not the same as for 
long time evolution. Short time instabilities are usually localized in space, while 
long time eigenvectors of (|23|l possesses a global structure. 

The third time scale corresponds to the error growing times more than 100 
days. On these scales the model exhibits either non- linear chaotic behavior with 
exponential growth of all X{t) (as it is the case in the square box), or stable behavior 
when a perturbation of initial state decreases with time (as it is the case in the Black 
sea model). 

In order to zoom these time scales, we plot the same data in the Log-Log and 
Log-Linear coordinates in figS] on the left and on the right respectively. One can 
see the error growth in the square box on this time scale is purely exponential with 
the same exponent A(t) = 74exp(0.027t). The multiplier A is particular for each 
parameter, but the exponent is always the same. This confirms the remark made 
in [13j, pT]: no matter how the perturbation was introduced into the model, it's 
long-time growth is determined by the model's dynamics. 

Comparing the evolution of the sensitivity of the model to different parameters, 
we see that on small scales it is the bottom topography that the model is the most 
sensitive to (thin solid line in figHj). An error in the topography produces 13 times 
bigger perturbation in the model state than a similar error in the model's initial 
conditions (thick solid line in figS]). However, A(t) does not grow at all on medium 
scales due to significant changes in the eigenvectors pattern. This leads to the fact 
that on long scales, the sensitivity of the model to the bottom topography is about 
2 times lower than the sensitivity to initial conditions. 

On the other hand, the sensitivity of the model to the discretization of operators 
near the boundary exhibits the opposite behavior. On short scales, corresponding A 
is 2 times lower than A obtained for perturbations of (po, but there is no stagnation of 
the growth on the middle scales. As a result, we see that the model is 4 times more 
sensitive to a than to (po for long error growing times. Moreover, small perturbation 
of initial conditions decreases in the Black sea configuration, while the perturbation 
of a results in an increasing perturbation of the solution. 

5 Modification of the boundary conditions 

As it has been noted, it is useless to analyze the set of obtained coefficients a to 
understand the modification of the boundary conditions. Instead of this, we shall 
see the difference between velocity fields with classical and with optimal coefficients 
a similarly to [21j. This difference has been averaged in time over 200 days time in- 
terval in order to reveal persistent modifications of the fiow produced by the optimal 
discretization. 

This average difference of the velocity together with the original velocity are 
presented in figiSJ We zoom the Southern part of the Black sea because it is in this 
region the difference shows the biggest values reaching 5 while in the middle of 
the sea it rarely exceeds 1 
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Figure 5: Original velocity field (left) and its modification (right) near the boundary 

We can note several principal features of the flow that have been modified by 
boundary conditions. First, we can see a strong current on the boundary. The slip 
condition (vanishing tangential velocity) has been replaced by a permanent current 
along the boundary. Moreover, impermeability condition has also been modified. 
The flow is now allowed to leave the domain ensuring, however, the global mass 
balance. One can see a strong persistent vortex centered at 42.2°N, 32.8°E which 
southern part crosses the boundary resulting in not only tangential but normal flux 
also. Similar vortices with lower amplitude can also be seen in places where the 
boundary changes direction. Optimal discretization allows the flow to cross the 
boundary in places where the direction change is not smooth. 

Tangential velocity component is amplified in the direct vicinity of the boundary. 
In these nodes we see a strong eastward flow that was forbidden by the boundary 
conditions in the classical formulation of the model. On the other hand, the eastward 
velocity is lower at nodes distanced by several grid cells from the boundary. At 
these nodes we see westward flow in the difference of the optimally discretized and 
classical models. That means the flow is moved towards the boundary, allowing 
more optimal representation of a thin current on a coarse grid that brings the model 
solution towards observations. 



6 Conclusion 

The comparative study presented in this paper shows the influence of different model 
parameters on the solution. The study is confined to the analysis of a low resolution 
model with a rather limited physics. Consequently we must acknowledge the results 
may be valid only in the described case. Additional physical processes (baroclinic 
dynamics, variable density due to heat and salinity fluxes, etc.) may modify results 
of this study revealing other parameters the model may be sensitive to. 

The main conclusion we can made from this comparison is the important role 
played by the boundary conditions on rigid boundaries. Almost all experiments 
show that the model is the most flexible with respect to control of coefficients a, 
this control allows us to bring the model's solution closer to the solution of the 
high-resolution model or to the observed data. 

Optimal a found in the assimilation window remain optimal long time after 
the end of assimilation improving the forecasting ability of the model. We could see 
that the fourth order model in the square box allows us to divide by two the forecast 
error of the 20 days forecast. Optimal a obtained in one month assimilation remains 
optimal even for a one year run of the Black sea model. 
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Finally, the long time sensitivity of the model's solution to a exceeds the sensi- 
tivity to almost all other parameters including the sensitivity to initial conditions. 
A perturbation of a of a given small norm results in a bigger perturbation of the 
model's solution than a perturbation of some other parameter of an equal norm. 

However, we could see that the influence of boundary conditions is only im- 
portant on long time scales, i.e. time scales that exceeds the characteristic time 
of the domain. In both experiments presented above the characteristic time was 
approximately equal to 5 days and in both experiments the sensitivity to a becomes 
important on scale longer than 5 days. On the other hand, on short scales, it is 
the bottom topography that influences the most the model's solution. Both in the 
Black sea and in the square box the sensitivity to topography is approximately 40 
times more important than the sensitivity to a. 

In addition to that, we should note that usually prescribed boundary conditions 
(impermeability and no-slip conditions have been used here as the initial guess for 
the minimization of the cost function) seem not to be optimal for the model. As 
we can see in figlSl modifying a we can bring the model much closer to the high 
resolution model or to the observational data. But, the numerical scheme is strongly 
modified in the assimilation process violating even impermeability condition. 

Taking into account an important influence of the numerical scheme that intro- 
duces boundary conditions into the model, it is reasonable to think about identi- 
fication of the optimal scheme by data assimilation process instead of prescribing 
classical boundary conditions. 

Acknowledgments. Author thanks Gennady Korotaev from Marine Hy- 
drophysical Institute, National Academy of Sciences of Ukraine for providing the 
model parameters and data for the upper layer model of the Black Sea. 
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